//----------------------------------*-C++-*----------------------------------//
/*!
 * \file   testDD_Equations.cc
 * \author Jeremy Roberts
 * \date   Thu Sep 8, 2011
 * \brief  Test of diamond difference equations.
 * \note   Copyright (C) 2011 Jeremy Roberts
 */
//---------------------------------------------------------------------------//
// $Rev:: 130                                           $:Rev of last commit
// $Author:: j.alyn.roberts@gmail.com                   $:Author of last commit
// $Date:: 2011-09-06 22:55:28 -0400 (Tue, 06 Sep 2011) $:Date of last commit
//---------------------------------------------------------------------------//

#include <iostream>
#include <iomanip>
#include <cstdio>
#include <cstdlib>
#include <cmath>

#include "test/Testing.hh"
#include "test/TestInput.hh"

#include "angle/Definitions.hh"
#include "angle/MomentsIndex.hh"
#include "angle/Moment_to_Discrete.hh"

#include "tran/Moments_Field.hh"
#include "tran/Psi_Cell_Field.hh"
#include "tran/CartMesh.hh"
#include "tran/Traits.hh"
#include "tran/Definitions.hh"

#include "../utilities/Constants.hh"
#include "../utilities/DBC.hh"
#include "../utilities/Soft_Equivalence.hh"
#include "../utilities/SP.hh"

// Equations
#include "tran/DD_1D_Equations.hh"
#include "tran/DD_2D_Equations.hh"
#include "tran/DD_3D_Equations.hh"

using namespace std;
using namespace slabtran;
using namespace util;
using namespace constants;

// 1-d tests.
int testDD_1D();
// 2-d tests
int testDD_2D();
// 3-d tests
int testDD_3D();

int main(int argc, char *argv[])
{
  Require(argc == 2);
  int test = atoi(argv[1]);
  cout << " ------ testDD_Equations ------ " << endl;
  int test_value = 0;
  int number_of_tests = 4;

  try
  {
    if (test == 1 || test == number_of_tests) test_value += testDD_1D();
    if (test == 2 || test == number_of_tests) test_value += testDD_2D();
    if (test == 3 || test == number_of_tests) test_value += testDD_3D();
  }
  catch (std::exception &err)
  {
    cout << "ERROR: While testing testDD_Equations, "
         << err.what()
         << endl;
  }
  catch ( ... )
  {
    cout << "ERROR: While testing testDD_Equations, "
         << "An UNKNOWN exception was thrown."
         << endl;
  }

  if (test > number_of_tests || test < 1)
    cout << " test " << test << " not available. " << endl;
  if (test_value)
    cout << " failed with value = " << test_value << endl;
  else
    cout << " passed." << endl;
  return test_value;
}

//===========================================================================//
// 1-D TESTS
//===========================================================================//
int testDD_1D()
{
  //----------------------------------------------------------------------//
  // SETUP FAKE PROBLEM
  //----------------------------------------------------------------------//

  // Make a simple 1-d mesh.
  Mesh<_1D>::SP_mesh mesh;
  mesh = new Mesh<_1D> (10.0, // width [cm]
                        8);   // meshes

  // Make a simple angular mesh.
  AngularMesh<_1D>::SP_angularmesh angularmesh;
  angularmesh = new AngularMesh<_1D> (4,         // quad order
                                      GAUSSLEG); // quad type

  // Make an M operator
  int legendre_order = 1;
  Moment_to_Discrete<_1D> M(legendre_order);
  M.build(angularmesh);

  // Make a dummy flux.
  Moments_Field<DD_1D_Equations, _1D>
    phi(mesh, legendre_order, 0.0);

  // Make a dummy psi.
  Psi_Cell_Field<DD_1D_Equations, _1D>
    psi(mesh, angularmesh->number_angles(), 0.0);

  // Make a dummy in and out psi.
  Mesh<_1D>::Space_Vector in_psi  = 1.2358; // just a scalar in 1-d
  Mesh<_1D>::Space_Vector out_psi = 0;      // just a scalar in 1-d

  // Dummy source for this mesh.
  double q = 1.0;

  // Dummy sigma for this mesh.
  double sigma = 1.0;

  // Use the first cell.
  int i = 0;


  //----------------------------------------------------------------------//
  // TEST 1-D DD EQUATIONS
  //----------------------------------------------------------------------//

  // Setup the equation.
  DD_1D_Equations eq(mesh);

  // We will go through each octant and ensure psi is computed
  // correctly.  Then, we see if psi is added correctly to the moments.

  //---------------------------------------------------------------------//
  // Check the first octant, ( 1, 1, 1).
  int o = 0;
  eq.set_octant((angularmesh->octant(o)).get_octant());

  // Use the first mu and weight in this octant.
  Octant<_1D>::Omega mu = (angularmesh->octant(o)).angles(0);
  double wt = (angularmesh->octant(o)).weight(0);
  int angle = angularmesh->index(o,0);
  TEST(angle == 0);

  // Generate reference values.
  double factor         = 2.0 * mu / (10.0/8.0);
  double psi_ref        = (q + factor * in_psi) / (sigma + factor);
  double out_psi_ref    = 2.0 * psi_ref - in_psi;
  //                      psi     * 0.5*(2*l + 1)*P_l(mu) * wt
  double phi0_ref       = psi_ref * 0.5*(2*0 + 1)*1.0     * wt;
  double phi1_ref       = psi_ref * 0.5*(2*1 + 1)*mu      * wt;

  // Solve.
  eq.solve(i, sigma, q, mu, wt, M, in_psi, phi, out_psi, psi, angle, true);

  // Check angular fluxes.
  TEST(soft_equiv(out_psi, out_psi_ref));
  TEST(soft_equiv(psi(i, 0, angle), psi_ref));
  TEST(soft_equiv(psi(i+1, 0, angle), 0.0));

  // Check moments.
  TEST(soft_equiv(phi(i, 0, 0), phi0_ref));
  TEST(soft_equiv(phi(i, 0, 1), phi1_ref));
  TEST(soft_equiv(phi(i+1, 0, 0), 0.0));

  //---------------------------------------------------------------------//
  // Check the second octant, ( -1, 1, 1).
  o = 1;
  eq.set_octant((angularmesh->octant(o)).get_octant());

  // Use the first mu and weight in this octant.
  mu = (angularmesh->octant(o)).angles(0);
  wt = (angularmesh->octant(o)).weight(0);
  angle = angularmesh->index(o, 0);
  TEST(angle == 2);

  // Update reference values.
  phi0_ref += psi_ref * 0.5*(2*0 + 1)*1.0 * wt;
  phi1_ref += psi_ref * 0.5*(2*1 + 1)*mu  * wt;

  // Solve.
  eq.solve(i, sigma, q, mu, wt, M, in_psi, phi, out_psi, psi, angle, true);

  // Check angular fluxes.
  TEST(soft_equiv(out_psi, out_psi_ref));
  TEST(soft_equiv(psi(i, 0, angle), psi_ref));
  TEST(soft_equiv(psi(i+1, 0, angle), 0.0));

  // Check moments.
  TEST(soft_equiv(phi(i,   0, 0), phi0_ref));
  TEST(soft_equiv(phi(i,   0, 1), phi1_ref));
  TEST(soft_equiv(phi(i+1, 0, 0), 0.0));

  return 0;
}

//===========================================================================//
// 2-D TESTS
//===========================================================================//
int testDD_2D()
{
  //----------------------------------------------------------------------//
  // SETUP FAKE PROBLEM
  //----------------------------------------------------------------------//

  // Make a simple 2-d mesh.
  Mesh<_2D>::SP_mesh mesh;
  mesh = new Mesh<_2D> (10.0, // width [cm]
                        8);   // meshes

  // Make a simple angular mesh.
  AngularMesh<_2D>::SP_angularmesh angularmesh;
  angularmesh = new AngularMesh<_2D> (4,         // quad order
                                      LEVELSYM); // quad type

  // Make an M operator
  int legendre_order = 1;
  Moment_to_Discrete<_2D> M(legendre_order);
  M.build(angularmesh);

  // Make a dummy flux.
  Moments_Field<DD_2D_Equations, _2D>
    phi(mesh, legendre_order, 0.0);

  // Make a dummy psi.
  Psi_Cell_Field<DD_2D_Equations, _2D>
    psi(mesh, angularmesh->number_angles(), 0.0);

  // Make a dummy in and out psi.
  Mesh<_2D>::Space_Vector in_psi(0.444, 0.555);
  Mesh<_2D>::Space_Vector out_psi(0.0, 0.0);

  // Dummy source for this mesh.
  double q = 1.0;

  // Dummy sigma for this mesh.
  double sigma = 1.0;

  // Use the following cell.
  int i = 3;
  int j = 5;
  int cell = 8 * 5 + 3; // 43
  TEST(mesh->convert(i, j) == cell);

  //----------------------------------------------------------------------//
  // TEST 2-D DD EQUATIONS
  //----------------------------------------------------------------------//

  // Setup the equation.
  DD_2D_Equations eq(mesh);

  //---------------------------------------------------------------------//
  // Check the first octant, ( 1, 1, 1).
  int o = 0;
  eq.set_octant((angularmesh->octant(o)).get_octant());

  // Use the first omega and weight in this octant.
  Octant<_2D>::Omega omega = (angularmesh->octant(o)).angles(0);
  double mu  = omega[MU];
  double eta = omega[ETA];
  double wt  = (angularmesh->octant(o)).weight(0);
  int angle  = angularmesh->index(o, 0);
  TEST(angle == 0);

  // Generate reference psi values.
  double factor_I    = 2.0 * mu  / (10.0/8.0);
  double factor_J    = 2.0 * eta / (10.0/8.0);
  double psi_ref     = (q + factor_I * in_psi[I]
                          + factor_J * in_psi[J]) /
                       (sigma + factor_I + factor_J);
  double out_psi_ref[] = {2.0 * psi_ref - in_psi[I],
                          2.0 * psi_ref - in_psi[J]};

  // Generate reference moment values.
  double phi0n0_ref  = psi_ref * inv_four_pi*(2*0 + 1)*1.0 * wt;
  double phi1m1_ref  = psi_ref * inv_four_pi*(2*1 + 1)*eta * wt;
  double phi1p1_ref  = psi_ref * inv_four_pi*(2*1 + 1)*mu  * wt;

  // Solve.
  eq.solve(i, j, sigma, q, omega, wt, M, in_psi, phi, out_psi, psi, angle, true);

  // Check angular fluxes.
  TEST(soft_equiv(out_psi[I], out_psi_ref[I]));
  TEST(soft_equiv(out_psi[J], out_psi_ref[J]));
  TEST(soft_equiv(psi(i,   j, 0, angle), psi_ref));
  TEST(soft_equiv(psi(i+1, j, 0, angle), 0.0));

  // Check moments.      c    u  l  m
  TEST(soft_equiv(phi(cell,   0, 0, 0), phi0n0_ref));
  TEST(soft_equiv(phi(cell,   0, 1,-1), phi1m1_ref));
  TEST(soft_equiv(phi(cell,   0, 1, 1), phi1p1_ref));
  TEST(soft_equiv(phi(cell+1, 0, 0, 0), 0.0));

  //---------------------------------------------------------------------//
  // Check the second octant, (-1, 1, 1).
  o = 1;
  eq.set_octant((angularmesh->octant(o)).get_octant());

  // Use the first omega and weight in this octant.
  omega = (angularmesh->octant(o)).angles(0);
  mu    = omega[MU];
  eta   = omega[ETA];
  wt    = (angularmesh->octant(o)).weight(0);
  angle = angularmesh->index(o, 0);
  TEST(angle == 3);

  // Update reference moment values.
  phi0n0_ref += psi_ref * inv_four_pi*(2*0 + 1)*1.0 * wt;
  phi1m1_ref += psi_ref * inv_four_pi*(2*1 + 1)*eta * wt;
  phi1p1_ref += psi_ref * inv_four_pi*(2*1 + 1)*mu  * wt;

  // Solve.
  eq.solve(i, j, sigma, q, omega, wt, M, in_psi, phi, out_psi, psi, angle, true);

  // Check angular fluxes.
  TEST(soft_equiv(out_psi[I], out_psi_ref[I]));
  TEST(soft_equiv(out_psi[J], out_psi_ref[J]));
  TEST(soft_equiv(psi(i,   j, 0, angle), psi_ref));
  TEST(soft_equiv(psi(i+1, j, 0, angle), 0.0));

  // Check moments.      c    u  l  m
  TEST(soft_equiv(phi(cell,   0, 0, 0), phi0n0_ref));
  TEST(soft_equiv(phi(cell,   0, 1,-1), phi1m1_ref));
  TEST(soft_equiv(phi(cell,   0, 1, 1), phi1p1_ref));
  TEST(soft_equiv(phi(cell+1, 0, 0, 0), 0.0));

  //---------------------------------------------------------------------//
  // Check the third octant, ( 1,-1, 1).
  o = 2;
  eq.set_octant((angularmesh->octant(o)).get_octant());

  // Use the first omega and weight in this octant.
  omega = (angularmesh->octant(o)).angles(0);
  mu    = omega[MU];
  eta   = omega[ETA];
  wt    = (angularmesh->octant(o)).weight(0);
  angle = angularmesh->index(o, 0);
  TEST(angle == 6);

  // Update reference moment values.
  phi0n0_ref += psi_ref * inv_four_pi*(2*0 + 1)*1.0 * wt;
  phi1m1_ref += psi_ref * inv_four_pi*(2*1 + 1)*eta * wt;
  phi1p1_ref += psi_ref * inv_four_pi*(2*1 + 1)*mu  * wt;

  // Solve.
  eq.solve(i, j, sigma, q, omega, wt, M, in_psi, phi, out_psi, psi, angle, true);

  // Check angular fluxes.
  TEST(soft_equiv(out_psi[I], out_psi_ref[I]));
  TEST(soft_equiv(out_psi[J], out_psi_ref[J]));
  TEST(soft_equiv(psi(i,   j, 0, angle), psi_ref));
  TEST(soft_equiv(psi(i+1, j, 0, angle), 0.0));

  // Check moments.      c    u  l  m
  TEST(soft_equiv(phi(cell,   0, 0, 0), phi0n0_ref));
  TEST(soft_equiv(phi(cell,   0, 1,-1), phi1m1_ref));
  TEST(soft_equiv(phi(cell,   0, 1, 1), phi1p1_ref));
  TEST(soft_equiv(phi(cell+1, 0, 0, 0), 0.0));

  //---------------------------------------------------------------------//
  // Check the fourth octant, (-1,-1, 1).
  o = 3;
  eq.set_octant((angularmesh->octant(o)).get_octant());

  // Use the first omega and weight in this octant.
  omega = (angularmesh->octant(o)).angles(0);
  mu    = omega[MU];
  eta   = omega[ETA];
  wt    = (angularmesh->octant(o)).weight(0);
  angle = angularmesh->index(o, 0);
  TEST(angle == 9);

  // Update reference moment values.
  phi0n0_ref += psi_ref * inv_four_pi*(2*0 + 1)*1.0 * wt;
  phi1m1_ref += psi_ref * inv_four_pi*(2*1 + 1)*eta * wt;
  phi1p1_ref += psi_ref * inv_four_pi*(2*1 + 1)*mu  * wt;

  // Solve.
  eq.solve(i, j, sigma, q, omega, wt, M, in_psi, phi, out_psi, psi, angle, true);

  // Check angular fluxes.
  TEST(soft_equiv(out_psi[I], out_psi_ref[I]));
  TEST(soft_equiv(out_psi[J], out_psi_ref[J]));
  TEST(soft_equiv(psi(i,   j, 0, angle), psi_ref));
  TEST(soft_equiv(psi(i+1, j, 0, angle), 0.0));

  // Check moments.      c    u  l  m
  TEST(soft_equiv(phi(cell,   0, 0, 0), phi0n0_ref));
  TEST(soft_equiv(phi(cell,   0, 1,-1), phi1m1_ref));
  TEST(soft_equiv(phi(cell,   0, 1, 1), phi1p1_ref));
  TEST(soft_equiv(phi(cell+1, 0, 0, 0), 0.0));

  return 0;
}

//===========================================================================//
// 3-D TESTS
//===========================================================================//
int testDD_3D()
{
  //----------------------------------------------------------------------//
  // SETUP FAKE PROBLEM
  //----------------------------------------------------------------------//

  // Make a simple 3-d mesh.
  Mesh<_3D>::SP_mesh mesh;
  mesh = new Mesh<_3D> (10.0, // width [cm]
                        8);   // meshes

  // Make a simple angular mesh.
  AngularMesh<_3D>::SP_angularmesh angularmesh;
  angularmesh = new AngularMesh<_3D> (4,         // quad order
                                      LEVELSYM); // quad type

  // Make an M operator
  int legendre_order = 1;
  Moment_to_Discrete<_3D> M(legendre_order);
  M.build(angularmesh);

  // Make a dummy flux.
  Moments_Field<DD_3D_Equations, _3D>
    phi(mesh, legendre_order, 0.0);

  // Make a dummy psi.
  Psi_Cell_Field<DD_3D_Equations, _3D>
    psi(mesh, angularmesh->number_angles(), 0.0);

  // Make a dummy in and out psi.
  Mesh<_3D>::Space_Vector in_psi(0.444, 0.555, 0.666);
  Mesh<_3D>::Space_Vector out_psi(0.0, 0.0, 0.0);

  // Dummy source for this mesh.
  double q = 1.0;

  // Dummy sigma for this mesh.
  double sigma = 1.0;

  // Use the following cell.
  int i = 3;
  int j = 5;
  int k = 5;
  int cell = 3 + 5 * 8 + 5 * 8 * 8; // 363
  TEST(mesh->convert(i, j, k) == cell);

  //----------------------------------------------------------------------//
  // TEST 3-D DD EQUATIONS
  //----------------------------------------------------------------------//

  // Setup the equation.
  DD_3D_Equations eq(mesh);

  //---------------------------------------------------------------------//
  // Check the first octant, ( 1, 1, 1).
  int o = 0;
  eq.set_octant((angularmesh->octant(o)).get_octant());

  // Use the first omega and weight in this octant.
  Octant<_3D>::Omega omega = (angularmesh->octant(o)).angles(0);
  double mu  = omega[MU];
  double eta = omega[ETA];
  double xi  = omega[XI];
  double wt  = (angularmesh->octant(o)).weight(0);
  int angle  = angularmesh->index(o, 0);
  TEST(angle == 0);

  // Generate reference psi values.
  double factor_I = 2.0 * mu  / (10.0/8.0);
  double factor_J = 2.0 * eta / (10.0/8.0);
  double factor_K = 2.0 * xi  / (10.0/8.0);
  double psi_ref  = (q + factor_I * in_psi[I]
                       + factor_J * in_psi[J]
                       + factor_K * in_psi[K]) /
                    (sigma + factor_I + factor_J + factor_K);
  double out_psi_ref[] = {2.0 * psi_ref - in_psi[I],
                          2.0 * psi_ref - in_psi[J],
                          2.0 * psi_ref - in_psi[K]};

  // Reference moment values.
  double phi0n0_ref = psi_ref * inv_four_pi*(2*0 + 1)*1.0 * wt;
  double phi1m1_ref = psi_ref * inv_four_pi*(2*1 + 1)*mu  * wt;
  double phi1n0_ref = psi_ref * inv_four_pi*(2*1 + 1)*xi  * wt;
  double phi1p1_ref = psi_ref * inv_four_pi*(2*1 + 1)*eta * wt;

  // Solve.
  eq.solve(i, j, k, sigma, q, omega, wt, M, in_psi, phi, out_psi, psi, angle, true);

  // Check angular fluxes.
  TEST(soft_equiv(out_psi[I], out_psi_ref[I]));
  TEST(soft_equiv(out_psi[J], out_psi_ref[J]));
  TEST(soft_equiv(out_psi[K], out_psi_ref[K]));
  TEST(soft_equiv(psi(i,   j, k, 0, angle), psi_ref));
  TEST(soft_equiv(psi(i+1, j, k, 0, angle), 0.0));

  // Check moments.      c    u  l  m
  TEST(soft_equiv(phi(cell,   0, 0, 0), phi0n0_ref));
  TEST(soft_equiv(phi(cell,   0, 1,-1), phi1m1_ref));
  TEST(soft_equiv(phi(cell,   0, 1, 0), phi1n0_ref));
  TEST(soft_equiv(phi(cell,   0, 1, 1), phi1p1_ref));
  TEST(soft_equiv(phi(cell+1, 0, 0, 0), 0.0));

  //---------------------------------------------------------------------//
  // Check the second octant, (-1, 1, 1).
  o = 1;
  eq.set_octant((angularmesh->octant(o)).get_octant());

  // Use the first omega and weight in this octant.
  omega = (angularmesh->octant(o)).angles(0);
  mu    = omega[MU];
  eta   = omega[ETA];
  xi    = omega[XI];
  wt    = (angularmesh->octant(o)).weight(0);
  angle = angularmesh->index(o, 0);
  TEST(angle == 3);

  // Update reference moment values.
  phi0n0_ref += psi_ref * inv_four_pi*(2*0 + 1)*1.0 * wt;
  phi1m1_ref += psi_ref * inv_four_pi*(2*1 + 1)*eta * wt;
  phi1n0_ref += psi_ref * inv_four_pi*(2*1 + 1)*xi  * wt;
  phi1p1_ref += psi_ref * inv_four_pi*(2*1 + 1)*mu  * wt;

  // Solve.
  eq.solve(i, j, k, sigma, q, omega, wt, M, in_psi, phi, out_psi, psi, angle, true);

  // Check angular fluxes.
  TEST(soft_equiv(out_psi[I], out_psi_ref[I]));
  TEST(soft_equiv(out_psi[J], out_psi_ref[J]));
  TEST(soft_equiv(out_psi[K], out_psi_ref[K]));
  TEST(soft_equiv(psi(i,   j, k, 0, angle), psi_ref));
  TEST(soft_equiv(psi(i+1, j, k, 0, angle), 0.0));

  // Check moments.      c    u  l  m
  TEST(soft_equiv(phi(cell,   0, 0, 0), phi0n0_ref));
  TEST(soft_equiv(phi(cell,   0, 1,-1), phi1m1_ref));
  TEST(soft_equiv(phi(cell,   0, 1, 0), phi1n0_ref));
  TEST(soft_equiv(phi(cell,   0, 1, 1), phi1p1_ref));
  TEST(soft_equiv(phi(cell+1, 0, 0, 0), 0.0));

  //---------------------------------------------------------------------//
  // Check the third octant, (-1, 1, 1).
  o = 2;
  eq.set_octant((angularmesh->octant(o)).get_octant());

  // Use the first omega and weight in this octant.
  omega = (angularmesh->octant(o)).angles(0);
  mu    = omega[MU];
  eta   = omega[ETA];
  xi    = omega[XI];
  wt    = (angularmesh->octant(o)).weight(0);
  angle = angularmesh->index(o, 0);
  TEST(angle == 6);

  // Update reference moment values.
  phi0n0_ref += psi_ref * inv_four_pi*(2*0 + 1)*1.0 * wt;
  phi1m1_ref += psi_ref * inv_four_pi*(2*1 + 1)*eta * wt;
  phi1n0_ref += psi_ref * inv_four_pi*(2*1 + 1)*xi  * wt;
  phi1p1_ref += psi_ref * inv_four_pi*(2*1 + 1)*mu  * wt;

  // Solve.
  eq.solve(i, j, k, sigma, q, omega, wt, M, in_psi, phi, out_psi, psi, angle, true);

  // Check angular fluxes.
  TEST(soft_equiv(out_psi[I], out_psi_ref[I]));
  TEST(soft_equiv(out_psi[J], out_psi_ref[J]));
  TEST(soft_equiv(out_psi[K], out_psi_ref[K]));
  TEST(soft_equiv(psi(i,   j, k, 0, angle), psi_ref));
  TEST(soft_equiv(psi(i+1, j, k, 0, angle), 0.0));

  // Check moments.      c    u  l  m
  TEST(soft_equiv(phi(cell,   0, 0, 0), phi0n0_ref));
  TEST(soft_equiv(phi(cell,   0, 1,-1), phi1m1_ref));
  TEST(soft_equiv(phi(cell,   0, 1, 0), phi1n0_ref));
  TEST(soft_equiv(phi(cell,   0, 1, 1), phi1p1_ref));
  TEST(soft_equiv(phi(cell+1, 0, 0, 0), 0.0));

  //---------------------------------------------------------------------//
  // Check the fourth octant, (-1,-1, 1).
  o = 3;
  eq.set_octant((angularmesh->octant(o)).get_octant());

  // Use the first omega and weight in this octant.
  omega = (angularmesh->octant(o)).angles(0);
  mu    = omega[MU];
  eta   = omega[ETA];
  xi    = omega[XI];
  wt    = (angularmesh->octant(o)).weight(0);
  angle = angularmesh->index(o, 0);
  TEST(angle == 9);

  // Update reference moment values.
  phi0n0_ref += psi_ref * inv_four_pi*(2*0 + 1)*1.0 * wt;
  phi1m1_ref += psi_ref * inv_four_pi*(2*1 + 1)*eta * wt;
  phi1n0_ref += psi_ref * inv_four_pi*(2*1 + 1)*xi  * wt;
  phi1p1_ref += psi_ref * inv_four_pi*(2*1 + 1)*mu  * wt;

  // Solve.
  eq.solve(i, j, k, sigma, q, omega, wt, M, in_psi, phi, out_psi, psi, angle, true);

  // Check angular fluxes.
  TEST(soft_equiv(out_psi[I], out_psi_ref[I]));
  TEST(soft_equiv(out_psi[J], out_psi_ref[J]));
  TEST(soft_equiv(out_psi[K], out_psi_ref[K]));
  TEST(soft_equiv(psi(i,   j, k, 0, angle), psi_ref));
  TEST(soft_equiv(psi(i+1, j, k, 0, angle), 0.0));

  // Check moments.      c    u  l  m
  TEST(soft_equiv(phi(cell,   0, 0, 0), phi0n0_ref));
  TEST(soft_equiv(phi(cell,   0, 1,-1), phi1m1_ref));
  TEST(soft_equiv(phi(cell,   0, 1, 0), phi1n0_ref));
  TEST(soft_equiv(phi(cell,   0, 1, 1), phi1p1_ref));
  TEST(soft_equiv(phi(cell+1, 0, 0, 0), 0.0));


  //---------------------------------------------------------------------//
  // Check the fifth octant, ( 1, 1,-1).
  o = 4;
  eq.set_octant((angularmesh->octant(o)).get_octant());

  // Use the first omega and weight in this octant.
  omega = (angularmesh->octant(o)).angles(0);
  mu    = omega[MU];
  eta   = omega[ETA];
  xi    = omega[XI];
  wt    = (angularmesh->octant(o)).weight(0);
  angle = angularmesh->index(o, 0);
  TEST(angle == 12);

  // Update reference moment values.
  phi0n0_ref += psi_ref * inv_four_pi*(2*0 + 1)*1.0 * wt;
  phi1m1_ref += psi_ref * inv_four_pi*(2*1 + 1)*eta * wt;
  phi1n0_ref += psi_ref * inv_four_pi*(2*1 + 1)*xi  * wt;
  phi1p1_ref += psi_ref * inv_four_pi*(2*1 + 1)*mu  * wt;

  // Solve.
  eq.solve(i, j, k, sigma, q, omega, wt, M, in_psi, phi, out_psi, psi, angle, true);

  // Check angular fluxes.
  TEST(soft_equiv(out_psi[I], out_psi_ref[I]));
  TEST(soft_equiv(out_psi[J], out_psi_ref[J]));
  TEST(soft_equiv(out_psi[K], out_psi_ref[K]));
  TEST(soft_equiv(psi(i,   j, k, 0, angle), psi_ref));
  TEST(soft_equiv(psi(i+1, j, k, 0, angle), 0.0));

  // Check moments.      c    u  l  m
  TEST(soft_equiv(phi(cell,   0, 0, 0), phi0n0_ref));
  TEST(soft_equiv(phi(cell,   0, 1,-1), phi1m1_ref));
  TEST(soft_equiv(phi(cell,   0, 1, 0), phi1n0_ref));
  TEST(soft_equiv(phi(cell,   0, 1, 1), phi1p1_ref));
  TEST(soft_equiv(phi(cell+1, 0, 0, 0), 0.0));

  //---------------------------------------------------------------------//
  // Check the sixth octant, (-1, 1,-1).
  o = 5;
  eq.set_octant((angularmesh->octant(o)).get_octant());

  // Use the first omega and weight in this octant.
  omega = (angularmesh->octant(o)).angles(0);
  mu    = omega[MU];
  eta   = omega[ETA];
  xi    = omega[XI];
  wt    = (angularmesh->octant(o)).weight(0);
  angle = angularmesh->index(o, 0);
  TEST(angle == 15);

  // Update reference moment values.
  phi0n0_ref += psi_ref * inv_four_pi*(2*0 + 1)*1.0 * wt;
  phi1m1_ref += psi_ref * inv_four_pi*(2*1 + 1)*eta   * wt;
  phi1n0_ref += psi_ref * inv_four_pi*(2*1 + 1)*xi  * wt;
  phi1p1_ref += psi_ref * inv_four_pi*(2*1 + 1)*mu  * wt;

  // Solve.
  eq.solve(i, j, k, sigma, q, omega, wt, M, in_psi, phi, out_psi, psi, angle, true);

  // Check angular fluxes.
  TEST(soft_equiv(out_psi[I], out_psi_ref[I]));
  TEST(soft_equiv(out_psi[J], out_psi_ref[J]));
  TEST(soft_equiv(out_psi[K], out_psi_ref[K]));
  TEST(soft_equiv(psi(i,   j, k, 0, angle), psi_ref));
  TEST(soft_equiv(psi(i+1, j, k, 0, angle), 0.0));

  // Check moments.      c    u  l  m
  TEST(soft_equiv(phi(cell,   0, 0, 0), phi0n0_ref));
  TEST(soft_equiv(phi(cell,   0, 1,-1), phi1m1_ref));
  TEST(soft_equiv(phi(cell,   0, 1, 0), phi1n0_ref));
  TEST(soft_equiv(phi(cell,   0, 1, 1), phi1p1_ref));
  TEST(soft_equiv(phi(cell+1, 0, 0, 0), 0.0));

  //---------------------------------------------------------------------//
  // Check the seventh octant, ( 1,-1,-1).
  o = 6;
  eq.set_octant((angularmesh->octant(o)).get_octant());

  // Use the first omega and weight in this octant.
  omega = (angularmesh->octant(o)).angles(0);
  mu    = omega[MU];
  eta   = omega[ETA];
  xi    = omega[XI];
  wt    = (angularmesh->octant(o)).weight(0);
  angle = angularmesh->index(o, 0);
  TEST(angle == 18);

  // Update reference moment values.
  phi0n0_ref += psi_ref * inv_four_pi*(2*0 + 1)*1.0 * wt;
  phi1m1_ref += psi_ref * inv_four_pi*(2*1 + 1)*eta * wt;
  phi1n0_ref += psi_ref * inv_four_pi*(2*1 + 1)*xi  * wt;
  phi1p1_ref += psi_ref * inv_four_pi*(2*1 + 1)*mu  * wt;

  // Solve.
  eq.solve(i, j, k, sigma, q, omega, wt, M, in_psi, phi, out_psi, psi, angle, true);

  // Check angular fluxes.
  TEST(soft_equiv(out_psi[I], out_psi_ref[I]));
  TEST(soft_equiv(out_psi[J], out_psi_ref[J]));
  TEST(soft_equiv(out_psi[K], out_psi_ref[K]));
  TEST(soft_equiv(psi(i,   j, k, 0, angle), psi_ref));
  TEST(soft_equiv(psi(i+1, j, k, 0, angle), 0.0));

  // Check moments.      c    u  l  m
  TEST(soft_equiv(phi(cell,   0, 0, 0), phi0n0_ref));
  TEST(soft_equiv(phi(cell,   0, 1,-1), phi1m1_ref));
  TEST(soft_equiv(phi(cell,   0, 1, 0), phi1n0_ref));
  TEST(soft_equiv(phi(cell,   0, 1, 1), phi1p1_ref));
  TEST(soft_equiv(phi(cell+1, 0, 0, 0), 0.0));

  //---------------------------------------------------------------------//
  // Check the eigth octant, (-1,-1,-1).
  o = 7;
  eq.set_octant((angularmesh->octant(o)).get_octant());

  // Use the first omega and weight in this octant.
  omega = (angularmesh->octant(o)).angles(0);
  mu    = omega[MU];
  eta   = omega[ETA];
  xi    = omega[XI];
  wt    = (angularmesh->octant(o)).weight(0);
  angle = angularmesh->index(o, 0);
  TEST(angle == 21);

  // Update reference moment values.
  phi0n0_ref += psi_ref * inv_four_pi*(2*0 + 1)*1.0 * wt;
  phi1m1_ref += psi_ref * inv_four_pi*(2*1 + 1)*eta * wt;
  phi1n0_ref += psi_ref * inv_four_pi*(2*1 + 1)*xi  * wt;
  phi1p1_ref += psi_ref * inv_four_pi*(2*1 + 1)*mu  * wt;

  // Solve.
  eq.solve(i, j, k, sigma, q, omega, wt, M, in_psi, phi, out_psi, psi, angle, true);

  // Check angular fluxes.
  TEST(soft_equiv(out_psi[I], out_psi_ref[I]));
  TEST(soft_equiv(out_psi[J], out_psi_ref[J]));
  TEST(soft_equiv(out_psi[K], out_psi_ref[K]));
  TEST(soft_equiv(psi(i,   j, k, 0, angle), psi_ref));
  TEST(soft_equiv(psi(i+1, j, k, 0, angle), 0.0));

  // Check moments.      c    u  l  m
  TEST(soft_equiv(phi(cell,   0, 0, 0), phi0n0_ref));
  TEST(soft_equiv(phi(cell,   0, 1,-1), phi1m1_ref));
  TEST(soft_equiv(phi(cell,   0, 1, 0), phi1n0_ref));
  TEST(soft_equiv(phi(cell,   0, 1, 1), phi1p1_ref));
  TEST(soft_equiv(phi(cell+1, 0, 0, 0), 0.0));

  return 0;
}
